% cleaned up and optimized update from extractTrajParams_v3_3


function r = extractTrajParams_v3_4(traj,stimwindow,fps,dosmooth,minimalProcess,integrationWindow)

%% set basic parameters first:
if nargin < 3 || isempty(fps)
    fps = 10;
end
if nargin < 4 || isempty(dosmooth)
    dosmooth = true;
end
if nargin < 5 || isempty(minimalProcess)
    minimalProcess = false;
end
doplot = false;
dilationthresh = 1;     % do not bother with dilation onset if magnitude of dilation is too small
dilationthreshZ = 3;    % do not bother with dilation onset if magnitude of the Z score is too small
onsetgofthresh = .65;

% reject if too many values are missing
if sum(isnan(traj))/length(traj) > (1/3)
    r.maxDilation = [nan,nan];
    r.stimonset = nan;
    r.stimonsetRange = nan(1,2);
    r.fwhm = nan(1,2);
    r.maxContractionVelocity = nan;
    r.totalDilation = nan;
    r.maxContraction = [nan,nan];
    r.peakTotalDilation = nan;
    r.peakFinalDilation = nan;
    return;
end

% generate extended window to look at other parameters; extension goes 4
% seconds after the end of the stimulus
if nargin < 6 || isempty(integrationWindow)
    extendedwindow = [stimwindow(1),round(stimwindow(2)+4*fps)];
else
    extendedwindow = integrationWindow;
end

if extendedwindow(2) > length(traj)
    extendedwindow(2) = length(traj);
end

% generate windows for the dilation in second second of the stimulus
% (around the peak) and the last second of the stimulus
stimwindow_initial = round(stimwindow(1)+fps*1):round(stimwindow(1)+fps*2);
stimwindow_later = round(stimwindow(2)-fps*1):round(stimwindow(2));

%% calculate derived trajectories:
baseline = traj(1:(stimwindow(1)-1));
trajZ = traj-nanmean(baseline);
trajZ = trajZ/nanstd(baseline);

if dosmooth
    trajsm = (smooth(traj,'sgolay',1))';            %better for onset time calculation
    trajsm2 = (smooth(traj,10,'rlowess'))';
    trajZsm = (smooth(trajZ,'sgolay',1))';
else
    trajsm = traj;
    trajsm2 = traj;
    trajZsm = trajZ;
end

%% calculate maximum dilation:
stimlength = (stimwindow(2)-stimwindow(1))/fps;
if stimlength > 3
    [maxdil,I] = max(trajsm2(stimwindow(1):stimwindow(2)));
    maxdilZ = max(trajZsm(stimwindow(1):stimwindow(2)));
else
    [maxdil,I] = max(trajsm2(extendedwindow(1):extendedwindow(2)));
    maxdilZ = max(trajZsm(extendedwindow(1):extendedwindow(2)));
end
I = I/fps;
r.maxDilation = [maxdil,I];

if maxdilZ < dilationthreshZ || maxdil < dilationthresh
    disp('subthreshold dilation')
    dilationSent = false;
else
    dilationSent = true;
end

%% calculate dilation onset:
% fit with sigmoid - this works if the dilation is pretty clean
if dilationSent && ~minimalProcess
    
    %x = 1:round((r.maxDilation(2)+.5+stimwindow(1)/fps)*fps);
    %x = (round(stimwindow(1)-2.5*fps)):round((r.maxDilation(2)+.5+stimwindow(1)/fps)*fps);
    %y = trajsm(x);
    x = (round(stimwindow(1)-1.5*fps)):round((r.maxDilation(2)+.5+stimwindow(1)/fps)*fps);
    x = x(x>0);
    x = x(x<length(traj));
    y = traj(x);
    %y = trajZsm(x);
    minOnsetFrame = stimwindow(1);
    [~, gof,onsetFrame] = fitStimOnsetSigmoid(x, y, minOnsetFrame,doplot);
    if gof.rsquare < onsetgofthresh
        % fit with a line:
        x = stimwindow(1):round((r.maxDilation(2)+.5+stimwindow(1)/fps)*fps);
        x = x(x<length(traj));
        y = trajsm(x);
        %y = trajZsm(x);
        onsetFrame2 = [fitStimOnsetLinear(x,y,doplot);nan;nan];
        r.stimonset = onsetFrame2;
    else
        r.stimonset = onsetFrame;
    end
    r.stimonset = (r.stimonset - stimwindow(1))/fps;
else
    r.stimonset = nan(3,1);
end

%% calculate integral of dilation
x = (extendedwindow(1):extendedwindow(2))/fps;
y = trajsm2(extendedwindow(1):extendedwindow(2));
r.totalDilation = trapz(x,y);

x = (stimwindow_initial(1):stimwindow_initial(2))/fps;
y = trajsm2(stimwindow_initial(1):stimwindow_initial(2));
r.peakTotalDilation = trapz(x,y);

x = (stimwindow_later(1):stimwindow_later(2))/fps;
y = trajsm2(stimwindow_later(1):stimwindow_later(2));
r.peakFinalDilation = trapz(x,y);

%% calculate contraction:
contractwindowstart = extendedwindow(1)+round(.5*fps);
[maxContraction,Icontraction] = min(trajsm2(contractwindowstart:extendedwindow(2)));
r.maxContraction = [maxContraction,Icontraction/fps];


function onsetFrame = fitStimOnsetLinear(x,y,doplot)
% this is a more robust fitting in noisier regimes:
temp = y-prctile(y,25);
i1 = find(temp>0);i1 = i1(1);
temp = y-prctile(y,75);
i2 = find(temp>0);i2 = i2(1);
x2 = x(i1(1):i2(1));
y2 = y(i1(1):i2(1));
p = polyfit(x2,y2,1);
onsetFrame = -p(2)/p(1);
if onsetFrame < x(1)
    onsetFrame = nan;       %set such that onset cannot be before stim onset
end
if doplot
    figure
    plot(x,y,'.')
    hold on
    x1 = linspace(x(1),x(end),100);
    y1 = polyval(p,x1);
    plot(x1,y1,'r')
    plot(onsetFrame,0,'ko')
    hold off
end
return;

function [fitresult, gof,onsetFrame] = fitStimOnsetSigmoid(x, y, minOnsetFrame,doplot)
% sigmoidal fit to the onset of the dilation - this is the most accurate IF
% the data is clean

[xData, yData] = prepareCurveData( x, y );

% Set up fittype and options.
ft = fittype( 'a + b/(1+exp(-(x-c)/d))', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [min(y), 0, minOnsetFrame, 0];
opts.StartPoint = [0.199983356523623 0.0420755183031353 0.634879508519415 0.691702153310096];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

thresholdDilationRelative = .1; % as a percent of max
threshDilation = (fitresult.b-fitresult.a)*thresholdDilationRelative + fitresult.a;
%threshDilation = (fitresult.b-fitresult.a)*.05 + fitresult.a;

onsetFrame = -log((1/thresholdDilationRelative)-1) - fitresult.d + fitresult.c;
%onsetFrame = fitresult.c - (log((fitresult.b/(threshDilation-fitresult.a))-1)*fitresult.d);
ci = confint(fitresult);

onsetFrame_lb = -log((1/thresholdDilationRelative)-1) - fitresult.d + ci(1,3);
onsetFrame_ub = -log((1/thresholdDilationRelative)-1) - fitresult.d + ci(2,3);
% onsetFrame_lb = ci(1,3) - (log((fitresult.b/(threshDilation-fitresult.a))-1)*fitresult.d);
% onsetFrame_ub = ci(2,3) - (log((fitresult.b/(threshDilation-fitresult.a))-1)*fitresult.d);
if onsetFrame < minOnsetFrame
    onsetFrame = nan(3,1);
else
    onsetFrame = [onsetFrame;onsetFrame_ub;onsetFrame_lb];
end
% Plot fit with data.
if doplot
    figure( 'Name', 'sigmoidal fit' );
    h = plot( fitresult, xData, yData );
    legend( h, 'y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
    % Label axes
    xlabel( 'x', 'Interpreter', 'none' );
    ylabel( 'y', 'Interpreter', 'none' );
    grid on
    hold on
    y = fitresult.a + fitresult.b./(1+exp(-(onsetFrame-fitresult.c)/fitresult.d));
    plot(onsetFrame,y,'ko')
end

return;





